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ABSTRACT 

Observations of chondritic meteorites and their ancestors, dust grains in protoplanetary disks, reveal 
the existence of strong, frequent heating events. One possible energy source for the heating that melts 
chondrules and anneals dust grains is the magnetic field that mediates the accretion flow, feeding off 
the vast reservoir of gravitational potential energy. In the absence of extremely spatially intermittent 
magnetic reconnection however, it has seemed unlikely that the dissipation of magnetic fields into heat 
in current sheets could reach the temperatures required to melt chondrules, T ^ 1800 K. In this paper, 
we show that there is hitherto unexplored dramatic behavior in protoplanetary disk current sheets 
triggered by the strong, positive relation between the temperature and the conductivity that can be 
understood as an electrical short. This is in opposition to the more commonly assumed resistivity 
increase in a magnetic reconnection region. The effect acts to focus the current sheets into even 
narrower, higher current and temperature regions. We lay out the basic principles of this behavior in 
this paper. 

Subject headings: Instabilities - Magnetic reconnection - Magnetohydrodynamics - Plasmas - Proto- 
planetary disks 



1. INTRODUCTION 

Observations of protostellar disks have revealed their 
integrated properties, including masses (Beckwith & Sar- 
gent 1993) and accretion rates (Hartmann 1998). Com- 
positional gradients in the dust (Van Boekel et al. 2004) 
can be detected, as well as the difference between pre- 
dominantly amorphous and crystalline mineral struc- 
tures at different radii (Waelkens et al. 1996; Malfait et 
al. 1998). Meanwhile, we have direct evidence of condi- 
tions in the protosolar disk. Laboratory measurements 
of textural, mineralogical, chemical, and isotopic proper- 
ties of meteoritic chondrules and calcium- aluminum rich 
inclusions (CAIs), as well as related high-temperature 
materials in comet samples (Brownlee et al. 2006; Zolen- 
sky et al. 2006; Nakamura et al. 2008; Simon et al. 2008), 
give strong constraints on their local formation history 
and environment. They represent melts condensed and 
cooled from temperatures of 1500-1800 K at rates of 
around 100-1000 K/hour (Radomsky and Hewins 1990; 
Lofgren and Lanier 1990; Connolly et al. 1998; Scott & 
Krot 2005; Ebel 2006): far faster than disk dynamical 
timescales, but far slower than the free-space cooling time 
of millimeter sized objects. 

The primary source of angular momentum transport 
in protoplanetary disks appears to be the conversion 
of orbital kinetic energy into magnetic energy through 
the magnetorotational instability (MRI; Velikhov 1959; 
Chandrasekhar 1961; Balbus & Hawley 1991, 1998), with 
gravitational instability playing a larger role at earlier 
times (e.g Lodato & Rice 2004). The ionization struc- 
ture of the disk midplane may introduce a magnetically 
dead zone (Gammie 1996) with reduced but non-zero tur- 
bulent viscosity (Fleming and Stone 2003; Oishi & Mac 
Low 2009), although its exact structure remains contro- 
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versial (Glassgold et al. 1997; Sano et al. 2000; Ilgner 
and Nelson 2006; Umebayashi & Nakano 2009; Turner & 
Drake 2009). 

MRI draws on the huge amount of energy contained in 
the differential rotation of the disk to drive magnetohy- 
drodynamical (MHD) turbulence. The turbulence will 
dissipate that energy into heat. However, the heating 
will occur intermittently, not uniformly. MHD turbu- 
lence forms current sheets (Parker 1972, 1994; Cowley et 
al. 1997) that dissipate energy at far greater than the 
average rate, and can provide locations for magnetic re- 
connection to occur. Romanova (2011) noted in pass- 
ing the volume filling nature of this mechanism (also see 
Romanova et al. 2011). Hirose & Turner (2011) used a 
moderate resolution simulation to demonstrate that such 
current sheets forming in the atmosphere above a dead 
zone can locally heat gas well above the radiative equi- 
librium temperature. 

Magnetic reconnection can occur in such current 
sheets, enhancing their local heating rate in small regions 
throughout them (For a review, see Yamada et al. 2010). 
Because local turbulent fiuctuations always introduce fi- 
nite perturbations to the symmetry of the system, only 
point reconnection events are expected (Boozer 2010). 
The speed of reconnection at these point events appears 
to be determined by Hall dynamics at the X-point rather 
than Ohmic resistivity (Pritchett 2001; Birn et al. 2001). 

However, measurements of reconnection events on the 
Sun have demonstrated that anomalous resistivities as 
much as 10^^ times even the Hall value must act to in- 
crease the rate of reconnection to the observed value (Lin 
et al. 2007). Lazarian and Vishniac (1999) suggested that 
this occurs in a turbulent reconnection layer. Recent 
three-dimensional simulations by Lapenta and Bettarini 
(2011) have used low dissipation, three-dimensional sim- 
ulations to conclusively demonstrate that such a volume- 
filling, turbulent reconnection layer naturally forms af- 
ter secondary tearing, kink, and Ray leigh- Taylor insta- 
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bilities completely destabilize the classical Sweet-Parker 
structure. Reconnection events such as these will locally 
heat the gas far beyond even the bulk values expected 
in the current sheets formed by the MRI, an idea also 
suggested in the disk context by King & Pringle (2010). 

This paper is primarily concerned with the behavior of 
current sheets in protoplanetary disks. As the local tem- 
perature climbs above ~ 900 K the dominant source of 
free electrons in a protoplanetary disk is thermal ioniza- 
tion of alkaline metals. The temperatures are still well 
below the ionization energy, so the argument of the ex- 
ponential in the Boltzmann term is quite small: thus the 
ionization fraction Xe = ^e/^n depends steeply on tem- 
perature. We have found that this results in startling 
new behavior, quite different from classical reconnection, 
with electrical short circuits effectively forming in these 
regions. This effect is the opposite of the more com- 
monly assumed anomalous resistivity that increases in 
the reconnection region (Krall & Liewer 1971; Sato & 
Hayashi 1979; Yamada et al. 2010). While our mecha- 
nism narrows current sheets, similarly to ambipolar dif- 
fusion (Brandenburg & Zweibel 1994), the mechanism is 
Ohmic resistivity rather than a drift of the charge carriers 
relative to the neutral gas. Our mechanism is also differ- 
ent from previous work on partially ionized reconnection 
because that focused on the transition to the collisionless 
regime (Malyshkin & Zweibel 2011; Zweibel et al. 2011). 

In this paper, we lay out the basic principles of this 
novel behavior, and explore the implications for heat- 
ing in dusty protoplanetary disks in more detail in a 
companion paper McNally et al. (2012, hereafter Paper 
II). Closely related behavior in planetary atmospheres 
was found by Menou (2012), who called it the thermo- 
resistive instability; though that term was also used by 
Price et al. (2012) for a system where the resistivity in- 
creases with temperature. 

In Sect. 2 we describe the physical principles at play. 
In Sect. 3 we lay out a numerical approach to modeling 
this behavior, whose results are described in Sect. 4, and 
discussed in Sect. 5. 

2. TEMPERATURE DEPENDENT RESISTIVITY 

2.1. Spatially- varying resistivity 

Consider the induction equation in the presence of 
Ohmic resistivity rj 



dB 

dt 



V x{UxB-r]J) . 



(1) 



where B is the magnetic flux, J the current, 77 the resis- 
tivity, and t time. If the resistivity is spatially uniform as 
normally assumed, then the resistive term can be rewrit- 
ten: 



SB 

— = V X (1/ X B) + rjV^B, 



(2) 



where the effect of the resistivity is to diffuse the mag- 
netic field. However, if the resistivity varies spatially, we 
must consider its spatial derivative: 



dB 

dt 



V x{U xB)^ r]V'^B - {Vr]) x J. (3) 



If the resistivity shows strong spatial variation, the final 
term in Equation (3) can dominate over the diffusive one. 



In the limit of a one- dimensional system, varying along 
X with the magnetic field pointing along we have 



dBy 

dt 



-S^ [va:By] + rjdlBy + da^rjd^By. (4) 



If the dxTjdxBy term is dominant over the diffusive term 
Tjd^By^ the Ohmic resistivity can act to steepen, rather 
than broaden, magnetic field gradients. 

This consideration of the spatial variation of the resis- 
tivity comes into full focus if the resistivity drops steeply 
with temperature, which can occur in a mostly neutral 
plasma in temperature ranges where one or more species 
is being thermally ionized. In this case, a local positive 
temperature perturbation that increases ionization will 
drive a positive current perturbation, which can in turn 
result in an increase in the local Ohmic heating. 

2.2. Steady State 

The base problem is one common to elementary electri- 
cal circuits: electrical shorts, albeit in the current driven 
regime. To put this on a quantitative basis, let us begin 
by considering a base steady-state one-dimensional sys- 
tem of length L, aligned with the x-axis and centered at 
X = 0, with magnetic structure given by 

(5) 
(6) 

0. (7) 



By\±L/2 — ^Bo/2 
Jz = dBy/dx 
dBy/dt = dx {r]Jz) 



Under the assumption of uniform rj = rjo we have simply 
Jz = Bq/L = Jq. Note that we have assumed that the 
velocity is Ux = 0. 

We now perturb the system, decreasing the resistivity 
to 77 = 77' = ^770 in a region of width AL = JL, centered 
at X = 0. Equations (5) - (7) then require 



(1 -S)Ji^ SJ[ = Jo 
Ji =r]J[, 



from which we can derive 



Ji = [5 + {l-d)fj] ' Jo, 



(8) 
(9) 

(10) 

where Ji is the value of Jz in the region with rj = tjq and 
J[ is the value of Jz in the region with rj = rj\ As long 
as 7j < 1, J{ > Jo, which is quite natural as current will 
preferentially flow in regions of low resistivity. 

More interesting is the Ohmic dissipation. The dissi- 
pation in the perturbed region 77 always exceeds that 
in the unperturbed region, as the electric fields in the two 
regions are equal by the requirement of a steady state. 
It also exceeds the dissipation in the base state 770^0 
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(11) 



This implies that even small decreases in the resistivity 
result in increased heating (effectively an electrical short) 
as long as 6 < 1/2. This size limitation is required to 
maintain sufficient residual current in the non-perturbed 
region that the electric field resistively generated there 
can be the effective voltage source for the short. 

In a partially ionized medium where the resistivity 
decreases with temperature because thermal ionization 
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produces increased charge carrier density, one would ex- 
pect the perturbed region to continue heating, reducing 
its resistivity further. Interestingly, this effect actually 
will reduce the total energy dissipated: the total cur- 
rent is constant, and we are reducing the resistivity of 
the medium the current flows through. While this effect 
generates local hot spots, it also decreases the total rate 
at which magnetic energy is dissipated into heat. 

2.3. Instability Analysis 

Section 2.2 gives a qualitative picture of a system that 
seems likely to experience an instability that would lead 
to an ever narrowing current sheet with increasing tem- 
perature and local heating. A real system will not be so 
idealized of course. We can explore the instability condi- 
tions quantitatively by performing a linear stability anal- 
ysis on a slightly less simplified one-dimensional model. 
We assume an incompressible fluid that cools (presum- 
ably radiatively) to a background bath temperature 
in a time t^. The equations are then 

-dx [-r]dxBy\ , 



dt 
dT 



TovidxByY 



(12) 
(13) 



where the Ohmic heating time is t^- The Ohmic heating 
term is normalized at the reference temperature Tq by the 
Ohmic heating with 77 = 770 and J = Jq. More exactly, 
we define the heating time 



th 



AnrinkBTo 
(7 -1)^70 Jo" 



(14) 



where the factor of Air comes from our use of cgs electro- 
magnetic units. Tin is the neutral number density and we 
assume a low ionization fraction. We track the temper- 
ature dependence of the resistivity through the general 
equation 



df] 
dT 



To 



(15) 



where Ti parameterizes the strength of the temperature 
gradient of 77. This also allows us to define the heating 
time of the resistivity tj^' through 



tnn 



To 


df] 




th 


df 


To 



Ti 
To 



th. 



(16) 



In a dominantly neutral medium in LTE, the Saha equa- 
tion is a simple approximation to the thermal ionization 
behavior, and the resistivity is dominated by the ioniza- 
tion fraction (see Equations 32 and 33 for more detail). 
In that case, considering only the exponential term in 
Equation 33 for analytic simplicity. Equation (15) be- 
comes 

dr] 
df 







To 


(-1) 



(17) 



where Ti is the temperature associated with first ioniza- 
tion . With Ti = 2.5188 x 10^ K, the ionization temper- 
ature of potassium, this approximates the ionization be- 
havior of protoplanetary disks at T 1000 K, as potas- 
sium has a low first ionization energy {ksTi) and suf- 
ficient abundances (Fromang et al. 2002). In this case. 



Equation (16) becomes 



To, 



tr]' — rpth- 



(18) 



With these approximations and definitions, we can de- 
rive the dispersion relation for a perturbation of the form 
^ikx+\t applied to a base state with 

By{x) = Jox (19) 
T{x) = To. (20) 
The linearized equations for the perturbations are 



XAB = 



TjoJotkAT 2a D 

— 770 rAE 

-Ll 



(21) 



AT 2 AT ToikAB 
XAT= + 2———. (22) 

tb trj' th Jo 

Solving Equations (21) and (22) we find 



1 _^ 1 _^ 1 

tr tb trp' 



1 



1 



tr, 



: 0, (23) 



where tr = k'^/r]o is the resistive time of the perturba- 
tion. As and tb are all positive, it follows that the 
perturbation can exhibit exponential growth (A > 0) if 
the constant term in Equation (23) is negative, i.e. if 



tr]' < tb, 



(24) 



which condition is independent of tr, unlike the usual 
situation for reconnection. While in the following sec- 
tion we will assume a resistivity profile that gives Equa- 
tion (17), the existence of the instability requires only a 
strong enough gradient of the temperature dependence 
of 77 (as measured through Ti and t^/). 

The independence of Equation (24) on tr arises in part 
because the resistivity plays an equal role in in the resis- 
tive time tr and the Ohmic heating time th and in part 
because the magnetic field transport is mediated purely 
through resistive effects in the imposed absence of veloc- 
ity. The steepening of magnetic field gradients through 
this instability and standard resistive spreading of mag- 
netic fields occur through the same resistivity operator. 
We emphasize this point: the transport of magnetic fields 
into the dissipation region is resistive in nature, rather 
than advective. Accordingly, this mechanism does not 
immediately struggle with the problem of exhaust that 
has bedeviled attempts to understand observed fast re- 
connection in the solar corona and elsewhere. 

The above analysis is a significant simplification of ac- 
tual physical systems even beyond its one- dimensional 
nature. In particular, we note that in a physical recon- 
nection region, the background current is not constant 
in time, and would be expected to decay resistively until 
the unstable modes have had time to grow. This oc- 
curs through the diffusive term in Equation (4) applied 
to the background current, which have been set to by 
our choice of background state, but which will not be 
negligible in general. 

As we will see, the assumption of a single physical 
length scale imposed in the above analysis also breaks 
down in practice, with the high temperature region nar- 
rowing over time in the non-linear regime. Further, 
MRI-active protoplanetary disks are compressible and 
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expected to have minimal plasma /3 ~ 1-10. The simula- 
tions described in the next section include these compli- 
cations by implementing terms such as the Lorentz force, 
which acts to compress a current sheet, and adiabatic 
heating and cooling. Finally, we note that our analysis 
of the current sheet has been done along the shortest di- 
mension, and that the current sheet will be much larger 
in the perpendicular, neglected dimensions. Because of 
this, the approximation of cooling to a bath temperature 
is a noteworthy oversimplification, and any cooling may 
act to expand the high temperature regions by radiative 
heating of the surroundings, as treated in more detail in 
Paper II. 

3. EQUATIONS AND NUMERICAL METHODS 

3.1. Numerical methods 

To model the dramatic behavior suggested by the lin- 
ear analysis, we have written a one dimensional code us- 
ing sixth order finite diff"erences on a logarithmic grid. 
We use implicit time integration with the CVODE pack- 
age (Hindmarsh et al. 2005). This allows us to follow the 
large, spatially limited variations in the resistivity. The 
logarithmically spaced grid allows us to push the bound- 
aries towards infinity while retaining resolution in the 
center of the current sheet. The number of grid points 
used in the simulations reported here varies from 500 
to 1000. Only the right half of the domain is included 
{x > 0) and a symmetrical inner boundary condition is 
used. 

3.2. Equations 

We solve the MHD fluid equations in 1.5 dimensions, 
including x-gradients and ^-components of vectors, with 
the one non-ideal term being the Ohmic resistivity. We 
use a somewhat more exact model of thermal ionization 
dominated by potassium. Although the linear analy- 
sis presented in Section 2.3 considered cooling, for this 
model we neglect cooling terms, deferring that additional 
complexity to Paper II. We use cgs units: magnetic field 
in Gauss and density in g/cm~^. 

With these approximations, the MHD equations be- 
come 



dp 
dt " 

dvx 
dt 

dBy 
dt 

dT 

'dt ' 



-dx (pvx) , 



(25) 



-VxdxVx dxP - - — dxBl + dxCsdxVx, (26) 

-dx [vxBy - r]{x)dxBy] , (27) 

-dx {Tvx) - ctPvx 

+ ^ {dxByf + ctCs {dxVxf , (28) 

where C,s is a shock viscosity included for stability. The 
shock viscosity is given by 



Cs Cs max(-S^v^Ax^)+ 



(29) 



where the constant Cs is taken as 10, max()+ denotes 
taking the maximum positive value over five grid points 
or zero otherwise, and Ax is the grid spacing. The equa- 
tion of state is that of an ideal gas and ct is the conver- 



sion factor between temperature and energy: 
P = nnksT, 

(7-1) 

Ct = . , 



(30) 
(31) 



where we use 7 = 7/5, is the neutral number density 
(assumed to dominate) and is the neutral molecu- 
lar mass. The resistivity associated with a dominantly 
neutral gas is given by Balbus & Terquem (2001) 



7] = 234T^/2x-^cmVs 



(32) 



and the ionization fraction = n^jrin^ under the as- 
sumption that the species being ionized is predominantly 
neutral and thermally ionized, becomes 



3/4 



X exp 



(33) 



where a is the fraction of the ionizing species to the total 
neutral population. In our canonical model, we consider 
only the thermal ionization of potassium (Fromang et al. 
2002). 

At the densities p ~ 10~^g/cm^, mean molecular mass 
[i = 2.33 and potassium fraction a = lO"*" of our canon- 
ical model, this equation breaks down at T > 1600 K 
when the potassium is significantly ionized. At higher 
temperatures in protoplanetary disks, other metals will 
also begin to contribute to the free electrons. At lower 
temperatures, below T ~ 1000 K, the ionization fraction 
from thermal processes becomes so low that in any as- 
trophysical system some non- thermal ionization source, 
such as ionizing stellar radiation or radionuclide decay, 
will dominate over thermal ionization. Even if they did 
not, the physical length scales required to achieve any 
MHD action in the presence of so high a resistivity be- 
come absurd. While for physical purposes. Equation (33) 
only applies in the temperature range 1000 K-1600 K, we 
will consider evolution at starting temperatures of 500 K 
and 2000 K to help test predictions about the strength 
of the gradient of the resistivity with respect to temper- 
ature. 

3.2.1. Initial and boundary conditions 

We consider initial conditions with Vx = 0^ po = 
10^ g cm~^, a = 10~^ and 

By{x) = Bq tanh(x/4), (34) 

which reproduces the magnetic field of a Harris (1962) 
current sheet. We denote the initial conditions at the 
box edge with the subscript 0, and use box widths large 
enough compared to £0 that and Bq are functionally 
identical. The density and temperature are then set to 
counterbalance the Lorentz force in the center assuming 
an adiabatic compression, is a control parameter that 
sets the total magnetic energy in the simulation. As we 
are interested in the ability of the magnetic field to heat 
the gas, instead of labeling runs with Bq^ we label them 
with the plasma beta: the ratio of thermal to magnetic 
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pressure 



= S7rpokBTo/{fimpBQ). 



(35) 



A value of /3o ~ 1 signifies an initial magnetic field energy 
that could raise the temperature throughout the box by 
~ 50% if converted directly to heat. The conversion will, 
however, be localized in our models. We list our control 
parameters in Table 1. 

We set the background temperature Tq, from which we 
derive the boundary gas pressure Pg^o = poksTo/ jj^rrip, 
and total pressure Pq = Pg^o + B'^/{^t^). We then derive 
the pressure, density and temperature profiles from the 
magnetic field profile given in Equation (34) 



Pg{x) =Po- B{xf/{8TT) 

T{x) = n{p,{x)/p,,,f-"'^ 



(36) 
(37) 
(38) 



This initial condition includes a resistivity minimum at 
the origin due to the increased temperature there. Note 
that Equation (33) is a decreasing function of the den- 
sity. If we were to use an isothermal hydrostatic initial 
condition, there would be an initial resistivity increase 
at the origin, which can split the instability into two, 
forming a swallowtail in a space-time diagram. In that 
regard, our adiabatic- hydrostatic initial condition is also 
gentler than a constant density-hydrostatic one due to 
the smaller spatial variation in the initial temperature. 

While in the analysis of Section 2.3 we assume a time- 
constant background temperature, the spatial variation 
of 77 may not dominate over the resistive diffusion of mag- 
netic field in Equation (4), especially at early time. This 
results in a decaying current density at the origin, until 
the instability has time to kick in. While we could ini- 
tialize our system with an inflow to confine the magnetic 
field, this would cause significant compressive heating. 

We make use of the symmetry in the problem along 
the mid-plane of the current sheet in order to solve only 
one half of the reconnect ion region. We use zero- gradient 
boundary conditions on the outer boundary (while push- 
ing them towards infinity) and symmetric/antisymmetric 

TABLE 1 

Run parameters 



^0 (K) 


(cm) 


B'^ (Gauss) 


^0 






3 


49.5 


500 


5 X 10^0 


3.4 


38.5 






4 


27.8 






5 


17.8 






5 


35.3 






5.5 


29.1 




5 X 10^ 


6 


24.5 


990 


7.1 


17.5 






10 


8.8 






15 


3.9 






20 


2.2 






7.5 


23.7 


1500 


2 X 10^ 


12.5 


8.5 






15 


5.9 






30 


1.5 






12.5 


11.4 


2000 


2.5 X 10^ 


19 


4.9 






25 


2.8 






50 


0.7 



boundary conditions as appropriate at the origin. 
4. RESULTS 

In Figure 1 we show the space-time evolution of the in- 
stability as a function of the ratio of thermal to magnetic 
pressure 

P = Sir pokTo/ifiiripB^), (39) 

with and without temperature dependent resistivity. In 
the third column we see the typical nature of the in- 
stability: a strong current sheet develops, shown by the 
narrowing of the magnetic field. When the temperature 
dependent ionization is turned off (column 4) the recon- 
nection region diffuses outwards normally. We also see 
a difference between this system and the idealized one 
of Section 2.3: the current density in the center spreads 
resistively at early times so the background current is 
not constant in time. The difference between columns 3 
and 4 shows the importance of treating the temperature 
dependence of the resistivity in this system. 

4.1. Growth Rate 

We can use the linear growth rate given by Equa- 
tion (23) to estimate the growth rate of the instability 
in the nonlinear regime reached in the simulations. In 
this regime, the value of the growth rate A(x,t) depends 
strongly on both position and time because of the varia- 
tion of the parameters, especially 77, but also the initial 
resistive spreading of the current density (see Figure 1). 
We extend to the nonlinear case by computing the value 
of A(0, t) at the center of the current sheet and examining 
it to see if it saturates to a constant value in the nonlin- 
ear regime that provides a reliable estimate of the actual 
growth rate. We define a timescale tc{t) = X~^{0,t) that 
we use to test this hypothesis. 

In the following computation of tc, we use i = 
msix{B) / J {x = 0) as the approximate, time- varying, 
width of the current sheet, taking the place of k~^ in 
the linear analysis, and 



tr,c 

t 



7]' ,C 



max(B2) ' 

" 7-1^' 
_ ^ 



(40) 
(41) 
(42) 

(43) 



where all spatially varying quantities are determined at 
X = 0. Solving Equation (23) for A using the above 
definitions, we find 



(lAc) 



1 



1 + 



1 + 



l2Ti 7 - 1 
T /3c 



/3c 



+ 



(44) 



T /3c 



1/2^ 



In Figure 2 we plot J(0), and tc> normalizing by the 
resistive time at the start of the simulation, = fr,c(i = 
0). We can see that tc{t) is reasonably well behaved 
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Fig. 1. — Space-time evolution of B in Gauss, self-normalized current density, and T in K for three different values of initial l3o, with a 
background temperature of To = 990 K. In the fixed rj case, the ionization level is set to that from the background non-thermal ionization. 
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Fig. 2. Top panels: /3o = 4.9; Middle panels: /5o = 18.6; 
Bottom panels: /3o = 37.5. Left panels: Time variation of J at 
the origin, normalized to its minimum value. Right panels: the 
dynamically defined timescale tc, normalized to the resistive time 
tr] . All for three values of /3o given in the panels and a background 
temperature of To = 990 K. 



Fig. 3. — Time series of J at the origin, normalized to its mini- 
mum value for two different time normalizations and six values of 
/3o- The background temperature of these models is To = 990 K. 
The solid, straight curves show ex.p(t/tc) and exp(t/2tc)- 



and indeed has a defined plateau that occurs after the 
onset of current density growth. On the other hand, 
it initially grows significantly because of the resistivity 
drop caused by Ohmic heating. As tc possesses a defined 
maximum value in systems that show current density 
growth, we will use its maximum value for our estimate 
of the growth rate. Unfortunately, this is not strictly well 
defined in cases that do not show current density growth 
(e.g. Figure 2, bottom panel). 

In Figures 3 and 4 we show the evolution of the cen- 
tral current density J for four temperatures and vary- 
ing f3o. Further, we include the functions exp(t/tc) and 
exp(t/2tc). If the analysis of Section 2.3 were exactly 
applicable with our definition of tc, then the growing in- 
stabilities would have the slope of the former. It is clear 
that, although tc is a reasonable estimate of the timescale 
for instability growth, the curves of growing J do not all 



have the same slope, and tc overestimates growth rates 
in some cases. Recall that tc is not well defined for the 
runs that fail to go unstable. 

The fact that tc is a good estimate of the growth rate 
of the instability is perhaps surprising in light of the be- 
havior of the resistivity (Figure 5, bottom panel). As 
the instability grows the resistivity at the origin drops 
by nearly two orders of magnitude, while the current 
sheet gets strongly concentrated: the system has become 
strongly nonlinear. A possible explanation for the con- 
tinued relevance of the linear analysis, however, is that 
the symmetry of the model problem maintains the va- 
lidity of the assumptions in Section 2.3. In particular, 
the current sheet concentrated by the instability has a 
flat current density at the origin that acts as the back- 
ground current density for further growth. We discuss 
this further in Section 4.3 below. 
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Fig. 4. — Time series of J at the origin, normalized to its min- 
imum value, for three background temperatures, Tq = 500, 1500, 
and 2000 K. 



4.2. Stability Criterion 

From Section 2.3 we might expect that all our simula- 
tions should be unstable due to the lack of cooling. While 
Figures 3 and 4 show runs that have not gone unstable, 
it is unclear whether this is due to the additional physics 
that we have added to the problem changing the condi- 
tion for instability, or merely inadequate run time (note 
that longer run times require larger boxes to exile the 
boundaries to infinity). The growth rate does appear to 
drop with increasing /3o even when normalized to t^ This 
slower growth may be due to the lower value of tr^d^c in 
the high /3o case (see Equation 44). Unlike in the linear 
analysis, the resistive time also acts to spread the back- 
ground current sheet and will decrease the instability's 
growth rate, and may halt it altogether. 

However, in Figure 3, top panel, the curve associated 
with /3o = 9.9 is just distinguishable from the left axis, 
while the curve associated with /3 = 4.9 is not. Clearly 
the increase in instability growth time with /3o is pro- 
nounced, with the lowest obviously unstable /3o = 25.6 
presented having its instability kick in at t ^ 5000t^. We 
expect that higher /3os, if unstable, would require even 
longer. Even if the higher /3o runs are eventually unsta- 
ble, it is unlikely to be a matter of practical concern in 



Fig. 5. — Self normalized current density, mass density, and resis- 
tivity, and temperature in Kelvins, for /3o = 3.9 and a background 
temperature of 990 K. Note that the density drops in the central 
region implying mass outflow from the central region even while 
magnetic flux is transported inwards (as traced by the growing 
current). 

a physical system. We attribute this to the longer inter- 
val in which the resistive spreading term in Equation (4) 
dominates over the instability term, resulting in an ever 
increasing and so an ever decreasing A. 

4.3. Saturation 

In the absence of cooling and ever decreasing resistiv- 
ity, it is not clear how the instability can saturate. So 
long as the current density J remains differentiable at 
the origin, symmetry requires that Q^^J = 0, so the ap- 
proximation assumed in Section 2.3 remains good. The 
saturated instability must maintain the constant back- 
ground current assumed. Further, as we can see in Fig- 
ure 5, while the system is compressible, it is not very 
compressible, with pT approximately constant. This sus- 
tained satisfaction of the linear instability criterion may 
explain why the time-scale estimate tc, which is deter- 
mined from a linear stability analysis, performs as well 
as it does despite the non-perturbative evolution of the 
system shown in Figure 5. 

In Figure 6 we show the evolution of the same system 
at two different resolutions, with the magnetic field plot- 
ted according to physical position on the left, and on the 
right, the current density plotted according to position on 
the logarithmically spaced grid. The upper panels show 
models with minimum grid resolution dx = 0.004 while 
the lower panels have dx = 0.002 ^o- The physical extent 
associated with the right panels is the same for the two 
resolutions and the same mapping of current density to 
color is used. In the left panels, we show the initial field, 
and the narrowing associated with the growing current 
sheet, which evolves similarly at both resolutions. In the 
right panels, we show the current density reaching the in- 
nermost grid point (the left side of the plot is the mirror 
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Fig. 6. — Magnetic field (Gauss), and current density plots for 
(3q = 2.2 and a background temperature of 990 K. Top panels: 
dx/io = 0.004; Bottom panels: dx/io = 0.002. On the right pan- 
els, the x-axes are in grid points rather than position (logarithmic 
spacing) but the plots cover the same physical extent. We can see 
the high resolution run continuing to narrow after the low resolu- 
tion run has hit the grid scale and started ringing, as shown by the 
sooty dashes at the very top of the top-right panel. 

image of the right) at low resolution, and associated ring- 
ing, while the high resolution run continues to narrow. 
It too will eventually narrow to below its grid resolution 
and be subject to the same numerical instability. 

In the absence of cooling, or changing physics, such 
as an alteration to the resistivity equation, it appears 
that the instability will not saturate. If the background 
current density is constant in time and different iable at 
the origin, symmetry requires that it have a constant 
component which is linearly unstable in the absence of 
cooling. Formally, the rf^'^B term in the induction 
equation (3) prohibits non-differentiable current densi- 
ties, so that the actual saturation mechanism must in- 
volve changing physics. 

One clear possibility is that the action of cooling in 
combination with a change in the temperature depen- 
dence of T] should halt the instability, as cooling will be- 
come stronger at higher temperatures while the resistiv- 
ity temperature gradient weakens. At high temperatures, 
it is expected that the temperature dependence of resis- 
tivity will change from that given by Equations (32) and 
(33). As the ionization fraction approaches full ioniza- 
tion, the temperature dependence of the resistivity will 
weaken. If the ionization fraction saturates, then the 
current sheet instability itself will saturate. 

5. DISCUSSION AND CONCLUSIONS 

We have shown that in slab-symmetric reconnection, a 
strong negative temperature dependence of the resistiv- 
ity can lead to an instability that concentrates the cur- 
rent in a narrow, high temperature, low resistivity sheet. 



This scenario is the polar opposite of the more common 
situation where the resistivity appears to increase inside 
current sheets through some anomalous resistivity (Krall 
& Liewer 1971; Sato & Hayashi 1979). Unlike many re- 
connection scenarios, the inward transport of the mag- 
netic field in our case is resistive rather than advective 
in nature, sidestepping issues of fluid pile up that oc- 
cur with advective field transport, as demonstrated by 
Figure 5, where the central density declines with time. 

However, rather than speeding the dissipation of mag- 
netic energy into heat, in one dimension the total dissipa- 
tion rate actually falls, thanks to the formation of small 
volumes with very low resistivity, even though heating 
increases sharply within those regions. This raises inter- 
esting questions about the structure of magnetic turbu- 
lence in three-dimensional systems with similar temper- 
ature dependent resistivities. 

We expect such systems to also show the concentra- 
tion of current into localized regions with high current 
and temperature (either two-dimensional sheets or one- 
dimensional tubes), with the rest of space taken up by 
almost force- free magnetic fields. As the concentrated 
current regions have low resistivity, they could poten- 
tially have long lifetimes, perhaps much longer than that 
associated with the high wavenumber tail of a subsonic 
turbulent cascade. Where the magnetic field energy does 
not exceed equipartition with the fluid kinetic energy, 
the bending of the magnetic field will create new cur- 
rent structures, leading to a highly balkanized magnetic 
field configuration, with the potential for extremely large 
and localized Lorentz forces because of the highly con- 
centrated current densities. 

While we have shown that, in sufficiently restricted 
circumstances, this instability occurs for any inverse re- 
lationship between the resistivity and the temperature, 
we have also shown that it can take a prohibitive time to 
set in. In practice it appears that the growth rate esti- 
mate of Equation (23), while sometimes an overestimate, 
is accurate within factors of a few. Unfortunately, evalu- 
ating it requires that the instability be growing, and any 
initial transients can result in strong overestimates of the 
growth rate (see Figure 2, early times). 

Although we have considered this effect from the per- 
spective of protoplanetary disks, it should occur in any 
system with an adequately strongly increasing conductiv- 
ity dependence on temperature when compared to avail- 
able cooling. Candidates include cool stellar surfaces 
with poorly ionized hydrogen, and even planetary atmo- 
spheres, as recently suggested by Menou (2012). 
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